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ABSTRACT 

Coalescing supermassive black hole binaries are produced by the mergers of galaxies and are the most powerful 
sources of gravitational waves accessible to space-based gravitational observatories. Some such mergers may 
occur in the presence of matter and magnetic fields and hence generate an electromagnetic counterpart. In this 
Letter, we present the first general relativistic simulations of magnetized plasma around merging supermassive 
black holes using the general relativistic magnetohydrodynamic code Whisky. By considering different mag- 
netic field strengths, going from non-magnetically dominated to magnetically dominated regimes, we explore 
how magnetic fields affect the dynamics of the plasma and the possible emission of electromagnetic signals. In 
particular we observe a total amplification of the magnetic field of ~ 2 orders of magnitude which is driven by 
the accretion onto the binary and that leads to much stronger electromagnetic signals, more than a factor of 10^ 
larger than comparable calculations done in the force-free regime where such amplifications are not possible. 
Keywords: accretion, accretion disks — black hole physics — gravitational waves — magnetohydrodynamic s 
(MHD) — methods: numerical 
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1. INTRODUCTION 

Space-based gravitational- wave (GW) detectors, such as 
the planned eLISA/NGO and SGO detectors, are expected 
to detect tens of supermassive black hole (BH) mergers per 
year. These detections will provide superbly precise measure- 
ments of the redshifted masses of the holes as well as the lu- 
minosity distance to the event. However, it is not possible 
to extract the redshift directly from the GWs. For this it is 
necessary to look for electromagnetic signatures that would 
identify the host galaxy. The resulting combination of the 
redshift with the luminosity distance would provi de a pow- 
erful cosmological probe ([ Hughes & Holz 2003; B erti et al. 
loops'; Ko csis et al.|2006t[Amn et al.|2009t .T: would also al- 
low precise tests of whether GWs travel at the speed of light, 
as required by general relativity. Although the merger itself 
produces no electromagnetic emission, if there are significant 
electromagnetic fields or mass near by in an ac cretion disk 
then there are various possibilities ( Schnittman| |201 1 ). For 
some disk accretion rates and binary mass ratios, the binary 
reaches a point in its coalescence such that further inspiral by 
emissio n of GWs occurs mo re rapidl y than the disk diffuses 
inward ([Armitage & Nataraj an 2002i |Milosavljevic & Phin-] 
ney |2U05p . This leads to a hole in the disk which is filled 
gradually after merger, leading to a source that b rightens over 
weeks to years depending on various parameters (Armitage &| 
Natarajan 2002; Milosavljevic & Phinney2005 ; Krolik 2010J 
Tanaka & Menou 2010; Shapiro 2010| ). Several authors have 



discussed consequences of the recoils from asymmetric emis 
sion of GWs during the coalescence, from pr ompt shocks to 



delayed emission lasting millions of year s ([Shields & Bon 
[ning 2008 ; Schnittman & Krolik 2008 ; Megevand et al.|2009t " 
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Lippai et al |200 ?, 'Corral es etal.|2010t [Anderson et al.|2010 
Rossi et al.|2010 ; Zanotti et al.|2010[ ) Emission might occui 



in the late inspiral from effects such as enhanced accretion, 
periodic Newtonian perturbation s, or shearing of the disk due 
to GWs ( [Kocsis & Loeb|[2QQ8 ). Earlier precursors are also 
possible, and in some cases the error volume from the GW 
signal may be small enough that the host galaxy can be iden- 
tified by morphology, mass, or by the presence of an active 
galactic nucleus. 

In the last few years there have been a number of publi- 
cations describing the evolution of gas and magnetic fields 
around merging supermassive BHs. van Meter and collab- 
orators performed test-particle simulations of the motion of 
accreting gas during the l ast phase of inspiral of comparable- 
mass supermassive BHs (van Meter et al. 2010). These sim- 
ulations suggested that near merger a significant fraction of 
particles can collide with each other at speeds approaching 
the speed of light, implying that a burst of radiation might ac- 
company the coalescence. Other works have instead started 
to investigate the effect that the merging BHs would have 
on surrounding gas and the possible emission of electromag- 
netic signals (O'Neill et al. 2009; Farris et al. 2010; Farris 
eraL||2011, ,Bode et al.,,2010, .Bogdanovic et al.,,2011, Bode 



et al. 2012). At the same tmie there have been the first m- 



vestigations of the effect of binary bla ck hole (BBH) mergers 
on electromagnetic fields i n vacuum ( [Palenzuela et al.|l2009| 
|2010c"; "Mosta et aL||2010| and in a niagne tically do minated 
pTasma (Pal enzuela et al.|2010b|al [Mosta etal.|2012| ). 
These studies have shown how magnetic and electric fields 



can be distorted by the motion of the BHs and hence lead 
to possible electromagnetic emiss ion. In particular, recent 
studies by P alenzuela and Mosta (Palenzuela e t al.[[2010b|a 



[Mosta et al. 2012) have raised the possibility that the motion 
of two BHs in a magnetically dominated plasma, i.e., in the 
so called force-free regime, could generate two separate jets, 
one around each BH, during the inspiral. At the time of the 
merger these two collimated jets would enter in contact and 
form a single jet emitted from the spinning BH formed af- 
ter the merger. It is, however, still unknown how general this 
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scenario is and whether the emission would be detectable. 

In this Letter, we present the first results from general 
relativistic magnetohydrodynamic (GRMHD) simulations of 
magnetized plasmas around merging supermassive black 
holes. By considering the evolution of equal-mass BBH sys- 
tems in plasmas with different levels of magnetization we fill 
the gap between the studies of non-magnetized gas and the re- 
sults obtained in the force-free and electro- vacuum regimes. 
We use a spacelike signature ( — ,+,+,+) and will typically 
use a system of units in which c = G = M = 1, where M 
is the total mass of the binary. In these units 1 M is equiv- 
alent to r^ 0.14 Mg hr and to -- 4.86 x 10"^ Mg pc, where 
Mg = M/(1O^M0). 

2. NUMERICAL METHODS AND INITIAL DATA 

Most of the details on the mathematical and numerical setup 
used for producing the results presented here are discussed in 
depth else where (Pollney et al.||2Q0 7t |Thomb urg||2Q0 4l p^ 
acomazzo & Rezzolla 2007; Giacomazzo et al. 2009, 2011 



Loffler et al.,2012j ). In what follows, we limit ourselves to a 
brief overview. 

2.1. Magnetohydrodynamics and Einstein Equations 

The evolution of the spacetime was obtained using the 
Ccatie code, a three-dimensional finite-differencing code 
providing the solution of a conformal tracele ss formulation 
of the Einstein equations ( [Pollney et al.||2007| ), and we used 
the ''mov ing puncture" method a nd gauge conditions devel- 
oped in van Mete r et al.| ( |2006| ). The GRMHD equations 
were instead solved using the Whisky co de (Giacomazzo, 
& Rezzolla| |2007t [Giacomazzo et al.|[2011|) , which adopts a 
vativ e form ulation of the equations as presented 
(|2006J) and high-resolution shock-capturing 
Its presented here have been computed 



flux-conservative 



Anton et al 

All the resu 



m 

schemes 

using the piecewise parabolic method, while the Harten-Lax- 
van Leer-Einfeldt appr oximate Riemann solver has been used 
to compute the fluxes ( Giacomazzo & Rezzolla|[2007 l. All 
the simulations were performed using a polytropic equation of 
state (EOS) with a polytropic exponent 7 = 4/3 and a poly- 
tropic constant n = {).2. We used a polytropic EOS instead of 
an ideal-fluid EOS because the computation of primitive from 
conservative variables in highly magnetized plasmas is much 
simpler and mo re robust for polytropes (see [Giacomazzo & 



|Rezzolla|2007| for details). The main results of our work are 
unaffected by this choice. 

In order to guarantee the divergence-free character of the 
M HD equations we evolve the vector potential as described 
in [Giacomazzo et al.|2011[ When evolving the vector poten- 
tial a gauge choi ce needs to be made and we here use the "al- 
gebraic gauge" ( Etienne et al.||2012 ) which was also us ed in 
previous GRMHD simulations with the Whisky code ( |Gia- 
[comazzo et al.|201 1 ; Rezzolla et al. 201 1 ). The code has been 
validated against a series of tests in special relativity ( Giaco 



mazzo & Rezzo llal2006l) and in full general relativity (Giaco- 
mazzo & Rezzolla 2007 ). 



Since the simulations performed here consider a plasma 
with a total mass negligible with respect to the mass of the 
two BHs, we have decoupled the Einstein equations from the 
matter dynamics, i.e., the metric variables are evolved using 
Einstein equations in vacuum. The same was done in t he gen- 
eral relativis tic hydrodynamic simulations reported in | Farris| 
et al.| ( 2QT0| ). Moreover, in order to prevent the formation of 
nonphysical values in the MHD quantities, we have excised 
the MHD variables inside the apparent horizon of each BH. 



2.2. Adaptive Mesh Refinement 

Both the Einstein and the GRMHD equations are solved 
using the vertex-centered adaptive mesh-refinement (AMR) 
approach provided by the Carpet driver (Schnetter et al. 
|2004| ). Our rather basic form of AMR consists of center- 
ing the highest-resolution level around each BH and in mov- 
ing the "boxes" following the position of the two BHs. For 
the results presented here we have used 11 refinement lev- 
els with the finest resolution being 0.0375M and the coarsest 
resolution being 38.4M. The finest grid has a radius of 6M 
whereas the coarsest grid extends to 1536 M. The large extent 
of our finest grid allows us to follow accurately the dynamics 
of the plasma around the BHs and it is also sufficiently large 
to avoid the spurious magnetic field amplifications that may 
occur when evolving the vector potential with the "algebraic 
gauge" (Etienne et al . 2012). 

For all the simulations reported here we have also used a 
reflection-symmetry condition across the z = plane and a 
TT-symmetry condition across the x = plane|^ 

2.3. Initial Data 

We have considered equal-mass systems of two non- 
spinning BHs with an initial separation of 8.48 M; these in- 
spiral for approximately three orbits before merger. The 
initial data were computed usi ng the public avai lable 
TwoPuncture code developed in |Ansorg et al.| ( [20Q4| ) and 
we chose the momentum of the punctures m order to en- 
sure that the orbit of the two BHs is quasicircular. We 
have considered two models: BO is surrounded by a non- 
magnetized plasma, and model B2 instead has an initially uni- 
form magnetic field with a ratio of magnetic to gas pressure 
(Pmag/Pgas = /^~^) initially equal to 2.5 x 10~ . The mag- 
netic field is aligned with the total angular momentum of the 
system, i.e., B^ = (0, 0, 5^) while in all the models the rest- 
mass density p is initially uniform and fills the entire domain; 
we have in mind a flo w that could be advection-dominated 
close to the holes (e.g., |Ichimaru|1977| and thus have a high 
enough radial velocity that it can keep up with the binary in- 
spiral throughout the entire coalescence. The initial distribu- 
tion of the magnetic field is similar to that adopted in previous 
works (Palenzue la et al. 2010b a , Mosta et a l. 2012) and it as- 
sumes that the magnetic field is anchored to a circumbinary 
disk located far outside of the numerical domain. We note 
that these initial data are similar t o what was used in previous 
non-magnetic (Farris et al."2010') and force-free (Palenzuela 
et al. 2010b a; Mosta et al. 2012) analyses. After this Letter 
was submitted. Noble et al. 2012 presented preliminary re- 
sults describing the evolution of a magnetized circumbinary 
accretion disk up to few orbits before the BH merger. Their 
simulations showed that magnetized plasma can indeed flow 
in the region where the BHs merge, and in future simulations 
we plan to use initial conditions more similar to the end re sult 
of the simulations recently reported in Noble et aL]p012] ). 



3. DYNAMICS 

Figure [T] shows the evolution of the rest-mass density p on 
the equatorial plane for models BO (left panels) and B2 (right 
panels). The evolution of the non-m agnetized model (lef t pan- 
els) resembles the one described in [Farris et aL] ( [201 0| ) with 
the formation of two spiral shocks during the inspiral and the 

^ Stated differently, we evolve only the region {x > 0, z > 0} and apply 
a 180° -rotational-symmetry boundary condition across the plane at x = 0. 
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Figure 1. Evolution of the rest-mass density p on the equatorial plane for the non-magnetized model BO (left panels) and for the magnetized case B2 (right 
panels). The units of time (shown at the bottom of each panel) and distance are M and the rest-mass density is normalized to its initial value. 
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Figure 2. Rest-mass density p on the xz plane for the non-magnetized model BO (left panel) and for the magnetized case B2 (right panel) at the end of the 
simulation {t ~ 470M). The unit of distance is M and the rest-mass density is normalized to its initial value. 
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Figure 3. Ratio of magnetic to gas pressure for the magnetized model B2 on the xy (left panel) and xz planes (right panel) at the end of the simulation 
{t ~ 470M). The unit of distance is M. Note that the color bar is different between these two panels, on the left panel any magnetically dominated region would 
be in red while in the right panel the minimum value of the color bar is Pmag/pgas = 1 in order to highlight only those regions that are magnetically dominated. 
The right panel contains five refinement levels; the finest refinement extends to z = 6 and the coarsest applies to z = 14 — 20. 
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Figure 4. Evolution of the Poynting vector on the xz plane for the magnetized model B2. The panels refer respectively to the time after two orbits, at the merger 
and at the end of the simulation. The units of time (shown at the bottom of each panel) and distance are M. 



formation of a central spinning BH surrounded by a spheri- 
cal distribution of matter accreting onto it. The magnetized 
model (right panels) shows quite different dynamics. Dur- 
ing the evolution the magnetic field strength increases by ap- 
proximately two orders of magnitude and contributes signif- 
icantly to the total pressure in the gas. Because of this jj the 
two shock waves that are present during the inspiral of model 
B are strongly reduced and hardly visible (first panel in the 
right column). Moreover, the density close to each of the two 
BHs and in the region connecting the two BHs is much larger 
than in the unmagnetized case. In this region it is also pos- 
sible to see the formation of instabilities that are not present 
without magnetic fields. Soon after the merger (second panel 
in the right column) the spinning BH is surrounded by a disk 
with a density a factor of ~ 3 larger than in the unmagnetized 
case, and two shock waves are formed and the system finally 
relaxes to its final configuration (third panel in the right col- 
umn). The temperature in the magnetic simulation is larger by 
up to ~ 40% for model B2 than in the unmagnetized model. 

Although the evolution on the equatorial plane already 
shows some differences due to the effect of the magnetic field, 
the main difference is in the rest-mass density on the merid- 
ional plane. In Figure [2] we show the rest-mass density p on 
the xz plane for models BO (left panel) and B2 (right panel) at 
the end of the simulation (t ~ 470M). Whereas the end result 
of the evolution of the unmagnetized model is a plasma accret- 
ing spherically onto a spinning BH, in the case of the magne- 
tized model B2 the system forms a thin accretion disk and a 
funnel is created around the spin axis of the BH. Although 
at this time no relativistic jet is emitted (the Lorentz factor is 
lower than ~ 3 at the end of the simulation), such emission 
might exist at later times (which are outside the scope of the 
present Letter). 

The difference between the rest-mass density distribution in 
the unmagnetized and magnetized cases can be better under- 
stood by looking at Figure [3] which shows the ratio of mag- 
netic to gas pressure (i.e., the inverse of the plasma parameter 
(3) for the magnetized model B2 on the xy (left panel) and xz 
planes (right panel) at the end of the simulation (t ^ 4 TOM). 
On the equatorial plane no region is magnetically dominated 
and the central region inside the disk has larger values of /3 
than in the initial conditions, but the xz plane shows clearly 
the presence of a strongly magnetically dominated region 



close to the spin axis of the BH with (3 ~ 10~^. During 
the inspiral and merger magnetic field lines are indeed com- 
pressed and twisted causing the magnetic field to be amplified 
of approximately two orders of magnitude. This highly mag- 
netized region is responsible for the creation of the thin disk 
and its funnel. If these simulations had been run for many 
orbital times (which would have required much larger com- 
putational resources), the amplification of the magnetic field 
would likely have been even larg er, due to effects li ke the 
magneto-rotational instability (Bal bus & Hawley|1991 ). 



4. ELECTROMAGNETIC EMISSION 

In Figure [4] we show the evolution of the z component of 
the Poynting vector on the xz plane for model B2. We show 
in particular its outgoing component after two orbits (first 
panel), at the time of the merger (second panel) and at the 
end of the simulation (third panel). One of the main differ- 
ences with respect to force-free simulations is that we do not 
observe the two strong and distinct jets originating from each 
BH and the "jet" that is emitted by the system propagates 
slowly into the medium surrounding the binary. Whereas 
in the force-free scenarios the jet propagates in a very low- 
density medium where the inertia of the plasma can be ne- 
glected, in our simulations the jet has to break through the 
infalling medium outside the binary. We also note that in our 
scenario the emission is mainly collimated and parallel to the 
angular momentum of the binary and to the spin of the final 
BH with no sign of the dominant quadrupolar emission that 
was observed in the recent force-free simulations of Mo sta eti 
al. (2012]). We also note that th e non-collimated emission dis- 
covered in |Mosta et al.| (2012) is larger than their collimated 
emission and it is ^ 600 times smaller than our luminosity. 

Finally, in Figure [5] we show the Poynting flux luminosity 
computed for moderB2 (blue solid line). We have rescaled 
our results to consider a binary system with a total mass 
of 10^ Mq and immersed in a plasma with a rest-mass den- 
sity p — 10~^^gcm~^. This corresponds to having an ini- 
tial magnetic field of ~ lO^G for model B2, w hich is also 
the same m agnetic field strength considered in |Palenzueia| 



^ See JMizuno et al.| i 2009t for an example of how the magnetic pressure 
can affect shock formation . 



et al. I (2010b). The luminosity is computed at a distance of 
z = lOM. Model B2 shows the characteristic increase in lu- 
minosity during the inspiral, with a peak corresponding to the 
time of the merger, followed by a drop-off of a factor ~ 2. 
This is qualitat ively similar to what is observed in force-free 
simulations ( |Palenzuela et al.|r2010b| ), but our actual lumi- 
nosities are considerably higher. This happens because in our 
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Figure 5. Evolution of the luminosity for the magnetized model B2 (blue 
solid line). The luminosity is computed at a distance z — lOM for a bi- 
nary system with a total mass of 10^ M©, an initial rest-mass density of 
10~iig cm""^, and an initial magnetic field of ~ lO^G. 

ideal GRMHD simulations the magnetic field is amplified of 
~ 2 orders of magnitude. So even when starting with a mag- 
netic field of ~ lO^G, the final configuration has a field of 
~ 10^ G. If we were to compare with force-free simulations 
starting w ith such a high field we would obtain similar lumi- 
nosities ( Palenzuela et al.||2Q10b| ). Indeed, we note that our 
luminosity at the end is ^ ly 'CTg s~ ^ and in the simulations 



reported in Palenzuela et al. (2010b) is ~ 10 '^ergs"^ for a 
magnetic field that is ~ 2 orders of magnitude lower than the 
one we have at the end of our simulations. Since the lumi- 
nosity Lem scales as 5^ our values are consistent with those 
one would observe in a force-free regime. This again high- 
lights one of the main differences between our simulations 
and those assuming a force-free regime since the beginning. 
Because of the accretion of the plasma onto the BHs, mag- 
netic field lines are compressed and twisted driving the large 
amplification we observe. In a force-free regime, the mag- 
netic field is decoupled from the dynamics of the matter and 
such large amplifications cannot be obtained. 

5. CONCLUSIONS 

We have presented the first numerical GRMHD simulations 
of magnetized plasmas around merging supermassive BBHs. 
We have for the first time investigated the role of magnetic 
fields in the plasma dynamics and filled the gap between the 
works that have considered non-magnetized gas and the re- 
sults obtained in the force-free and electro- vacuum regimes. 

We have shown that even plasmas that are initially not 
magnetically dominated have different dynamics than in the 
unmagnetized case and that magnetic plasmas can generate 
strong and collimated electromagnetic emission. We therefore 
generalize the physical regimes of matter and electromagnetic 
fields around coalescing BHs that can lead to potentially de- 
tectable emission. 



We thank Phil Armitage, Tamara Bogdanovic, Bernard 
Kelly, Krzysztof Nalewajko, Carlos Palenzuela, Luciano Rez- 
zolla, Jeremy Schnittman, and Roman Shcherbakov for use- 
ful comments and suggestions. We also thank Philip Cow- 



perthwaite for help in visualizing some of the numerical data. 
Resources supporting this work were provided by the NASA 
High-End Computing (HEC) program through the NASA Ad- 
vanced Supercomputing (NAS) Division at Ames Research 
Center and NASA Center for Climate Simulation (NCCS) at 
Goddard Space Flight Center. Numerical simulations were 
also performed on the cluster RANGER at the Texas Ad- 
vanced Computing Center (TACC) at The University of Texas 
at Austin through XSEDE grant no. TG-PHYl 10027. B.G. 
acknowledges support from NASA grant no. NNX09AI75G 
and NSF grant no. AST 1009396. J.B. and J.v.M. acknowl- 
edge support from NASA grant 09-ATP09-0136. 

REFERENCES 

Anderson, M., Lehner, L., Megevand, M., & Neilsen, D. 2010, Phys. Rev. D, 

81,044004 
Ansorg, M., Brugmann, B., & Tichy, W. 2004, Phys. Rev. D, 70, 064011 
Anton, L., Zanotti, O., Miralles, J. A., et al. 2006, ApJ, 637, 296 
Armitage, R J., & Natarajan, R 2002, ApJ, 567, L9 
Arun, K. G., Mishra, C. K., Van Den Broeck, C, et al. 2009, Class. 

Quantum Grav., 26, 094021 
Balbus, S. A., & Hawley, J. R 1991, ApJ, 376, 214 

Berti, E., Buonanno, A., & Will, C. M. 2005, Class. Quantum Grav., 22, 943 
Bode, T., Bogdanovic, T., Haas, R., et al. 2012, ApJ, 744, 45 
Bode, T., Haas, R., Bogdanovic, T., Laguna, P., & Shoemaker, D. 2010, ApJ, 

715, 1117 
Bogdanovic, T., Bode, T., Haas, R., Laguna, P., & Shoemaker, D. 2011, 

Class. Quantum Grav., 28, 094020 
Corrales, L. R., Haiman, Z., & MacFadyen, A. 2010, MNRAS, 404, 947 
Etienne, Z. B., Paschalidis, V., Liu, Y. T., & Shapiro, S. L. 2012, 

Phys. Rev. D, 85, 024013 
Farris, B. D., Liu, Y. T., & Shapiro, S. L. 2010, Phys. Rev. D, 81, 084008 
Farris, B. D., Liu, Y. T., & Shapiro, S. L. 2011, Phys. Rev. D, 84, 024024 
Giacomazzo, B., & RezzoUa, L. 2006, J. Fluid Mech., 562, 223 
Giacomazzo, B., & Rezzolla, L. 2007, Class. Quantum Grav., 24, S235 
Giacomazzo, B., Rezzolla, L., & Baiotti, L. 2009, MNRAS, 399, L164 
Giacomazzo, B., Rezzolla, L., & Baiotti, L. 2011, Phys. Rev. D, 83, 044014 
Hughes, S. A., & Holz, D. E. 2003, Class. Quantum Grav., 20, 65 
Ichimaru, S. 1977, ApJ, 214, 840 

Kocsis, B., Frei, Z., Haiman, Z., & Menou, K. 2006, ApJ, 637, 27 
Kocsis, B., & Loeb, A. 2008, Phys. Rev. Lett., 101, 041101 
Krolik, J. H. 2010, ApJ, 709, 774 
Lippai, Z., Frei, Z., & Haiman, Z. 2008, ApJ, 676, L5 
Loffler, F., Faber, J., Bentivegna, E., et al. 2012, Class. Quantum Grav., 29, 

115001 
Megevand, M., Anderson, M., Frank, J., et al. 2009, Phys. Rev. D, 80, 

024012 
Milosavljevic, M., & Phinney, E. S. 2005, ApJ, 622, L93 
Mizuno, Y, Zhang, B., Giacomazzo, B., et al. 2009, ApJ, 690, L47 
Mosta, P., Alic, D., Rezzolla, L., Zanotti, O., & Palenzuela, C. 2012, ApJ, 

749, L32 
Mosta, R, Palenzuela, C, Rezzolla, L., et al. 2010, Ph ys. Rev. D, 81, 064 017 
Noble, S. C, Mundim, B. C, Nakano, H., et al. 2012, a rXiv: 1204. 1073| 
O'Neill, S. M., Miller, M. C, Bogdanovic, T., Reynolds, C. S., & 

Schnittman, J. D. 2009, ApJ, 700, 859 
Palenzuela, C, Anderson, M., Lehner, L., Liebling, S. L., & Neilsen, D. 

2009, Phys. Rev. Lett., 103, 081101 
Palenzuela, C, Garrett, T., Lehner, L., & Liebling, S. L. 2010a, 

Phys. Rev. D, 82, 044045 
Palenzuela, C, Lehner, L., & Liebling, S. L. 2010b, Science, 329, 927 
Palenzuela, C, Lehner, L., & Yoshida, S. 2010c, Phys. Rev. D, 81, 084007 
Pollney, D., Reisswig, C, Rezzolla, L., et al. 2007, Phys. Rev. D, 76, 124002 
Rezzolla, L., Giacomazzo, B., Baiotti, L., et al. 2011, ApJ, 732, L6 
Rossi, E. M., Lodato, G., Armitage, R J., Pringle, J. E., & King, A. R. 2010, 

MNRAS, 401, 2021 
Schnetter, E., Hawley, S. H., & Hawke, I. 2004, Class. Quantum Grav., 21, 

1465 
Schnittman, J. D. 2011, Class. Quantum Grav., 28, 094021 
Schnittman, J. D., & Krolik, J. H. 2008, ApJ, 684, 835 
Shapiro, S. L. 2010, Phys. Rev. D, 81, 024019 
Shields, G. A., & Bonning, E. W. 2008, ApJ, 682, 758 
Tanaka, T., & Menou, K. 2010, ApJ, 714, 404 
Thornburg, J. 2004, Class. Quantum Grav., 21, 743 



GRMHD SIMULATIONS OF MERGING SUPERMASSIVE BLACK HOLES 



van Meter, J. R., Baker, J. G., Koppitz, M., & Choi, D.-I. 2006, 

Phys. Rev. D, 73, 124011 
van Meter, J. R., Wise, J. H., Miller, M. C, et al. 2010, ApJ, 711, L89 
Zanotti, O., Rezzolla, L., Del Zanna, L., & Palenzuela, C. 2010, A&A, 523, 

A8 



